Direct volume rendering with high visual quality has now become a standard precondition for the visualization of three-dimensional medical data sets. The principle of volume rendering is already well known in the prior art. Volume rendering is used to determine a two-dimensional projection, here the display data set, of a three-dimensional source data set, including a display parameter, discretely sampled in voxels. In the medical field, this may entail anatomical parameters (e.g. Hounsfield values in a CT data set). However, other variables, in particular scalar variables, are also conceivable (e.g. variables determined from functional imaging and the like).
Direct volume rendering is based on the premise that each value of the display parameter is projected onto an extinction density and a color density. This act is called classification and uses a transfer function that may be formulated in two parts as an extinction function that projects a value of the display parameter onto an extinction density or opacity density, and a color function that projects a value of the display parameter onto a color density. Since the display parameter is present in discrete voxels, interpolation may be performed before or after the classification in order to obtain the smoothest possible curve over the volume acquired in the three-dimensional source data set. Therefore, either the densities or the display parameter may be interpolated.
Volume rendering processes may be based on a design eye point and a direction of view so that a ray results. The integration of an integrand containing the extinction function and the color function over the ray produces an intensity value describing the reproduction of the corresponding pixel in the display data set. The integral described is known as the volume-rendering integral and, from a mathematical viewpoint, has to be resolved exactly or at least approximately to determine the display data set, either explicitly or implicitly.
In its simplest form, the volume-rendering integral is resolved approximately by discretization using a Riemann sum, wherein the volume of the source data set is sampled along rays and the source data set is used as an input data record in order to determine color density values and extinction density values for the rays. If the sampling distance, and therefore the step size, of the Riemann sum converges toward 0, the discrete approximation also converges toward the exact solution of the volume-rendering integral.
Particularly in medical applications, if possible, the display data sets may be available in real time. Real time availability of the display data sets is difficult when so-called rendering servers are used in order to provide display data sets for multiple users. Therefore, the sampling distances used for approximations of this kind with Riemann sums are rather large. When using the Riemann sum, the resultant direct volume rendering images may display significant discretization artifacts. Although the number of these artifacts drops when the sampling distance is reduced, the sampling distance required to prevent such artifacts is dependent on how quickly the values of the display parameter in the volume change locally and on how quickly the transfer function changes over the value range of the display parameters in the voxels.
Therefore, the sampling distance required to avoid local artifacts is dependent on both the spatial frequencies of changes in the source data set and on the spatial frequencies of the transfer function so that the overall sampling distance required may be understood as a product of these two contributory parts. Therefore, an extremely small sampling distance is required in order to obtain artifact-free display data sets as a resultant images so that the performance of the rendering process on the computing mechanism carrying it out is extremely low. In addition, it is extremely difficult to predict the optimum sampling distance so that no artifacts occur.
To resolve this problem, a pre-integration technique was suggested. Standard publications suggesting pre-integration include the articles by N. Max et al., “Area and volume coherence for efficient visualization of 3d scalar functions”, ACM Computer Graphics (Proceedings San Diego Workshop on Volume Visualization 1990) 24, 5 (1990), pages 27-33 (hereinafter Max), S. Röttger et al., “Hardware-accelerated volume and isosurface rendering based on cell-projection”, Proceedings IEEE Visualization 2000 (2000), pages 109-116 (hereinafter Röttger1) and K. Engel et al., “High quality pre-integrated volume rendering using hardware accelerated pixel shading”, Proceedings Graphics Hardware 2001 (2001), Mark W., Schilling A., (Eds.), ACM Press, pages 9-16 (hereinafter Engel).
Pre-integration enables much higher image quality without having to increase the sampling rate and therefore reduce the sampling distance. A pre-integration table is generated as an n-dimensional table that uses the transfer function, in advance. Therefore, the extinction function and the color function are used as input values. For a linear change between all possible front and back positions of an interval, the volume-rendering integral is solved exactly and stored in the table. During the actual process of volume rendering, the data in the pre-integration table is used instead of the original transfer function because the volume-rendering integral may be broken down into contributions that may be read from the pre-integration table. This approach produces a marked increase in image quality compared to other approaches that calculate the Riemann sum explicitly with larger sampling distances.
Consideration of a one-dimensional transfer function results in a three-dimensional pre-integration table with two dimensions relating to the front and back values of the display parameter for an interval. The integration is performed over the sampling distance. Sampling distance dictates the third dimension storing the step size. However, the need to produce a high-quality three-dimensional pre-integration table for 16-bit volume data places large requirements on the required storage space. Moreover, it is important for a user to be able to change the transfer function interactively. For an high-quality, three-dimensional pre-integration table, interactive changing of the transfer function requires an immense amount of computing capacity and an immense amount of memory space in order to generate the pre-integration table in the desired time.
In cases that the sampling distance is assumed to be constant, the pre-integration table is only two-dimensional. Two-dimensional pre-integration tables are considered to be a practically applicable variable. However, assuming a constant sampling distance may not be possible many specific applications, even for simple cases where volume rendering is performed for an orthographic projection (“orthographic camera”).
The difficulty of assuming a constant sampling distance is demonstrated by the example of generating a display data set parts of a ray arriving at or leaving the volume covered by the three-dimensional source data set. Although it would be possible using a constant sampling distance or sampling rate to avoid one of these situations in that the sampling position is moved correspondingly along the ray, it is not conceivable for both boundary points to be treated correctly. Other examples include volume clipping and volume segmentation; both are simple extensions of volume rendering that are frequently used in post-processing applications for the evaluation of the medical data. Without the adaptation of the sampling distance or sampling rate, clearly identifiable rendering artifacts would occur in such situations.
The option of pre-integrated volume rendering with a discretionary sampling distance may be advisable when the actual three-dimensional source data sets are taken into consideration. For example, when trilinear interpolation is performed, the resultant volume is only C∞-constant within a group of eight, trilinearly interpolated voxels, but not on the boundary of the eight voxels where there is only C0-constancy. It would now be particularly advantageous for the sampling positions to be positioned at such C0-boundaries if the interpolated values are approximated along a ray. An approach frequently used for the optimization of performance is adaptive volume rendering. With adaptive volume rendering, the sampling distance may vary in dependence upon local properties.
The ability to use pre-integration with discretionary sampling distances is also extremely advantageous when gradient opacity modulation is to be performed with activated pre-integration. Gradient opacity modulation is a rendering technique enabling the modification of opacity using the local gradient strength during the rendering. Gradient opacity modulation is often used as an additional volume-rendering submode, enabling certain properties to be emphasized in dependence on the gradient strength, e.g. in order to emphasize border regions of tissues.
While the pre-integration technique represents a basic precondition for high-quality direct volume rendering, as yet no optimum algorithm is available in order to implement discretionary sampling distances. Although a few approaches have already been suggested that will be described below, these existing approaches either construe a three-dimensional pre-integration table, entail compromises with respect to the image quality of the display data set, or require a large number of accesses to pre-integration tables that is also disadvantageous with regard to performance on the computing mechanism.
Before a more detailed description of the prior art relating to freely selectable sampling distances with pre-integration, a brief introduction to the basic mathematical formulae for pre-integrated volume rendering is provided to enable better understanding relating to this subject. More detailed information on this subject is made in, for example, the article by M. Kraus and T. Ertl, “Pre-integrated volume rendering”, The Visualization Handbook (2004), Hansen C. D. Johnson C. R., (Eds.), Academic Press, pages 211-228 (hereinafter Kraus1).
Let τ(s) be a scalar function describing the extinction density as a function of the display parameter, s, e.g. the extinction function, and c(s)τ(s) a scalar function describing a color density (c(s) may also be termed a color function), the volume-rendering integral along a ray parameterized by z in the interval [0,D] is obtained as Equation 1:I=∫0Dτ(s(z)(c(s(z))exp(−∫0zτ(s(t))dt)dz  Eq. 1
where s(z) is a function that assigns an interpolated voxel value of the display parameter to the parameter z.
This well-known integral equation represents the mathematical basis for volume rendering. A volume-rendering implemented as hardware and/or software on a computing mechanism will evaluate this formula as some kind of approximation.
If the interval along the ray is divided into n equidistant parts with a sampling distance d=D/n, the opacity αi and the at least one color coefficient Ci for sub-intervals of the sampling distance, d, may be defined as Equation 2 and Equation 3:αi=1−exp(−∫id(i+1)dτ(s(t))dt)  Eq. 2Ci=∫id(i+1)dτ(s(z))c(s(z))exp(−∫idzτ(s(t))dt)dz  Eq. 3
When these sub-integrals are used, volume-rendering Equation 1 becomes Equation 4:I=Σi=0n−1CiΠj=0i−1(1−αj)  Eq. 4
This breakdown that is also the basis for a pre-integration technique, wherein in particular the Ci and the αi, the latter also optionally as the transmission Ti=1−αi, may provide the pre-integration tables and the same applies for the embodiments described below. Therefore, Equation 4 discloses the breakdown of the volume-rendering integral discussed above.
If, in one interval of the sampling distance, d, the function s(z) is approximated as linear, wherein, the corresponding relationship is determined by, in addition to the sampling distance, d, the front value, sf, and the back value, sb, of the display parameter specified by s(z), Equations 2 and 3 may be approximated as Equation 5 and 6:αi≈1−exp(−∫01τ((1−ω)sf+ωsb)ddω  Eq. 5Ci≈∫01τ((1−ω)sf+ωsb)c((1−ω)sf+ωsb)exp(−∫0ωτ((1−t)sf+tsb)ddt)ddω  Eq. 6
The problem-free transition to the integration variable ω enabled by the assumption of linearity results in an interval distance in Equations 5 and 6 that is no longer d but 1.
Equation 5 and Equation 6 are now only dependent on three parameters, namely the front value, sf, and back value, sb, of a sub-interval and on the sampling distance, d, of the sub-interval. If a three-dimensional pre-integration table is generated for all possible constellations of α(sf, sb, d) and C(sf, sb, d), the volume-rendering integral may be approximated by Equation 4, retrieving the corresponding opacities and color coefficients from the pre-integration table. Therefore, there is a breakdown into the sub-intervals that defines d, the values for sf and sb are determined at the sampling positions and the pre-integration table may be used.
Pre-integration may be considered an important contribution in volume rendering because resultant image quality of the display data set for a given sampling rate may be significantly increased compared to other approaches while the algorithm remains simple and quick as long as the pre-integration table is available. Without pre-integration, the high sampling rate required in order to prevent rendering artifacts would quickly result in performance problems for non-linear transfer functions. Accordingly, pre-integrated volume rendering used by virtually all volume-rendering providing high quality in computing mechanisms that are able to provide new display data sets in real time with interactive frame rates.
The publication on pre-integration by Röttger1, suggests an approach generating a three-dimensional pre-integration table and using three-dimensional texture hardware in order to save the pre-integration table. The pre-integration table is thus made available for inquiries. However, in many cases, a three-dimensional table of this kind may not be used in practice since high requirements are placed on the corresponding storage mechanism, and a long time or a high computing capacity is required in order to update the table. Therefore, the prior art has already suggested approaches that directly or indirectly avoid three-dimensional pre-integration tables and nevertheless may enable freely selectable sampling distances.
One category of approaches avoiding the use of three-dimensional pre-integration tables is the restriction of transfer functions to cases where the resultant integrals of the volume-rendering equation may be treated analytically. For example, in “Gaussian transfer functions for multi-field volume visualizations”, Proceedings Visualization 2003 (2003), pages 497-504, (hereinafter Kniss). Kniss suggested special transfer functions based on Gaussian functions. However, since approaches of this kind are not possible for general transfer functions, they are of less interest.
Engel suggests that Equation 5 and Equation 6 are approximated using one-dimensional integral functions. For the approximation of Equation 6, they suggested the use of an integral function that ignores the self-attenuation term resulting in a much simpler integrand, but this is based on the assumption that the expression τ(s)*d is very small. Another approximation that uses integral functions but makes fewer assumptions was suggested in a publication by M. Kraus, “Pre-integrated Volume Rendering for Multi-Dimensional Transfer Functions”, IEEE/EG Symposium on Volume and Point-Based Graphics, H.-C. Hege, D. Laidlaw, R. Pajarola, O. Staat, (Eds.), pages 1-8, 2008 (hereinafter Kraus2). In this publication, self-attenuation is not ignored but is also taken into account. This enabled the approximation to be reduced for general transfer functions and greater sampling distances.
The approach suggested by Engel, that the volume-rendering equation be simplified with the aid of integral functions, was originally suggested in order to accelerate the calculation of two-dimensional pre-integration tables and was later extended by Kraus in order to approximate pre-integrations for multidimensional transfer functions as well. While, therefore, the original purpose was not to deal with the case of discretionary sampling distances, these approaches are still suitable for the approximation of this case as well since the sampling distance, d, appears as an independent parameter in their equations that approximate the color coefficients using integral functions.
Problematically, the use of such approaches with integral functions is subject to limitations. For example, it is permissible to assume that the extinction density may be infinitely large such that the ray reaches full opacity directly when a point of this type is reached. In such a case, a one-dimensional integral function would not work well since the integral function is no longer able to reproduce opacity changes that occur after such a point of infinite opacity. Then, approaches with integral functions may only be used to a restricted degree if extremely high quality approximations are required.
In the following, let αi and Ci be the opacities and color coefficients for a sampling distance, d, and αi′ and Ci′ the opacities and color coefficients for the same front and back values but relating to a different sampling distance, d′. A category of algorithms is known use a two-dimensional pre-integration table providing C(sf, sb, d) and αi(sf, sb, d) for the sampling distance, d. The values C′ and α′ are to be determined for another discretionary sampling distance d′ but still from the two-dimensional pre-integration table. Examples algorithms are provided in the publications by S. Röttger et al., “Smart Hardware-Accelerated Volume Rendering”, Proceedings of EG/IEEE TCVG Symposium on Visualization VisSym '03 (2003), pages 231-238 (hereinafter Röttger2), and J. P. Schulze et al., Integrating pre-integration into the shear-warp algorithm”, Proceedings Volume Graphics 2003 (2003), pages 109-118 (hereinafter Schulze).
It is useful that the calculation of the opacity coefficient αi′ may be performed without an approximation, because the relationship expressed in Equation 7:1−αi′=(1−αi)d′/d  Eq. 7
is applicable between αi and αi′ (see e.g., Röttger2). Unfortunately, there is no corresponding relationship for the color coefficients. A coarse estimation ignoring self-attenuation was disclosed in Röttger1 by Equation 8:
                              C          i          ′                ≈                              C            i                    ⁢                                    ⅆ              ′                        ⅆ                                              Eq        .                                  ⁢        8            
If Equation 8 is used to calculate color coefficients for other sampling distances, d′, strong artifacts occur in conjunction with a darkening of the colors in the display data set. For this reason, a mathematical equation was derived in Schulze that is more suitable for approximating the correct color coefficients. The Schulze equation considers the special case of a volume-rendering integral with constant extinction density and the color density. If the solution of this special case is used as an approximation of the general case, the following approximation is obtained as Equation 9:
                                          C            i            ′                    ≈                                    C              i                        ⁢                                          α                i                ′                                            α                i                                                    =                              C            i                    ⁢                                    1              -                                                (                                      1                    -                                          α                      i                                                        )                                                                      ⅆ                    ′                                    ⁢                                      /                    ⅆ                                                                                      1              -                              (                                  1                  -                                      α                    i                                                  )                                                                        Eq        .                                  ⁢        9            
for the case when αi is unequal to zero and the equality does not apply for general color and extinction functions. If αi is equal to zero, the application of L'Hospital's rule to Equation 9 again results in Equation 8.
The advantageous feature of this approach is that the general pre-integration case for discretionary sampling distances is approximated by an algorithm that only requires a two-dimensional pre-integration table instead of a three-dimensional pre-integration table. However, this simple embodiment has drawbacks. The use of this approach entails the occurrence of a wide variety of rendering artifacts, in particular in regions where the opacity is high. This is due to the fact that this approach is only an approximation for the case of general color and extinction functions.
There is an interesting relationship between the approaches that use an integral function and the correction Equation 8 and Equation 9. The approximations for the color coefficients that are based on the integral functions, may also be used to derive equations permitting the calculation of opacities and color coefficients for a second sampling distance from opacities and color coefficients for a first sampling distance. The use of the integral formulae of Engel, obtains the color coefficients for Equation 8, while Krauss2 results in Equation 9.
A publication by Guthe et al., “High-quality unstructured volume rendering on the PC platform”, Proceedings Graphics Hardware 2002 (2002), pages 119-125 (hereinafter Guthe), disclosed an algorithm permitting high-quality color correction for general sampling distances. This algorithm utilizes the fact that one single pre-integration table is not sufficient to achieve high-quality correction of the color coefficients for different sampling distances. Guthe suggests the calculation of several pre-integration tables for normalized chromaticity Ci/αi for different fixed sampling distances and the approximation of the case of general sampling distances d′ by the use of an interpolation polynomial and scaling with the corrected opacity αi. Thus, the approximation of the color coefficients for the sampling distance d′ requires the generation of n two-dimensional pre-integration tables for normalized chromaticity for different sampling distances and the retrieval of values from all these n tables in order to interpolate them with a polynomial of the degree (n−1). The result is then scaled with the opacity α′ that, according to the other approaches, may be calculated directly by using Equation 7.
According to Guthe, the minimum n requiring this approach in order to achieve a good approximation quality is at least n=3 that is, at least three pre-integration tables have to be generated and, with volume rendering, it is necessary to make retrievals from all the tables and perform interpolation. A further increase of n further increases the mathematical quality of the numerical approximation, wherein, however, there is a simultaneous increase in the number of pre-integration tables required and the number of accesses to the tables.
Since according to Guthe, the normalized chromaticity is interpolated and the result corrected with the corrected opacity α′, this approach may be considered to be an extension of the approach according to Schulze with the latter being understood as the trivial case n=1. In this case, only the normalized chromaticity for the sampling distance d would be corrected with the corrected opacity, α′, resulting directly in Equation 9.
To summarize, at the present time at least three color integration tables are required in order to determine high-quality color coefficients for the case of general sampling distances, placing high requirements on storage space and access possibilities for the tables. A large degree of effort is required to determine the tables and provides further room for improvement with respect to quality.